#clean up
rm(list=ls())

#load libraries
library(akima)
library(geoR)
library(scatterplot3d)  
library(ggmap)
library(lattice)
library(foreign)
library(maps)
library(maptools)
library(spdep)
library(rgdal)
library(prevR)
library(car)
library(xtable)
library(rgdal)
library(rgeos)

#options
options(scipen=8)

#Set to working directory. Adjust the command below to reflect YOUR working directory.
#setwd('/Users/jamie/Documents/BIGdata/censusBlocks/')

###ALL DATA ARE FROM http://www2.census.gov/geo/tiger/TIGER2010BLKPOPHU/###
#Accessed 12 December 2015

###Load data in a loop###
codes<-c('01','02','04','05','06','08','09','10','11','12','13','15','16','17','18','19','20','21','22','23','24','25','26','27','28','29','30','31','32','33','34','35','36','37','38','39','40','41','42','44','45','46','47','48','49','50','51','53','54','55','56')

geoBlock<-matrix(NA,nrow=1,ncol=7)
colnames(geoBlock)<-c('STATEA','COUNTYA','TRACTA','BLOCKA','blockRadius','blockEast','blockNorth')

for(i in 1:51){
#for(i in 1:2){
	print(i)
	#load and project map
	file<-paste("tabblock2010_",codes[i],"_pophu/tabblock2010_",codes[i],"_pophu.shp",sep="")
	polys.0<-readShapePoly(fn=file,proj4string=CRS("+proj=longlat +datum=NAD83"))
	rm(file)
	polys.2<-spTransform(polys.0,CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=km +no_defs"))
	polys.3<-as(polys.2,'SpatialPolygons')

	#save radius and easting/northing centroids
	area<-gArea(polys.3,byid=TRUE)
	radius<-sqrt(area/pi)
	cent<-coordinates(polys.2)
	temporary<-cbind(polys.0@data[,1:4],radius,cent)
	colnames(temporary)<-c('STATEA','COUNTYA','TRACTA','BLOCKA','blockRadius','blockEast','blockNorth')
	geoBlock<-rbind(geoBlock,temporary)
	}
#head(temporary)
#head(geoBlock)

#write out data
geoBlock<-geoBlock[-1,]
write.csv(geoBlock,"geoBlock.csv",row.names=F)
